Empirical oscillating potentials for alloys from ab-initio fits 
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oscillating pair potentials" (EOPP) have the form 

C C 
V(r) = — + — C0S (k*r + cfi*) 
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By fitting to a database of ab-initio forces and energies, we can extract pair potentials for alloys, with a simple 
six-parameter analytic form including Friedel oscillations, which give a remarkably faithful account of many 
complex intermetallic compounds. As examples we show results for (crystal or quasicrystal) structure prediction 
and phonon spectrum for three systems: Fe-B, Al-Mg-Zn, and Al-Cu-Fe. 
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Several kinds of problem in materials modeling can be ad- 
dressed only by classical interatomic potentials. Consider 
complex alloy crystal structures: it is impossible to use 
diffraction-data refined structures straightforwardly, for they 
almost always contain sites with mixed or fractional occu- 
pancies. These must be resolved properly to evaluate phys- 
ical properties, which can be unfeasible with current fast ab- 
initio codes using density functional theory, such as VASP [1]. 
Only moderate resources are needed for single evaluations of 
the total energy, but repeated evaluations are prohibitive for 
unit cells having up to 10 3 atoms per cell; even if just 1-5% 
of the sites are uncertain, one must examine a vast number 
of variant structures to assign the occupancies optimally. In 
any case, classical potentials are required when the structure 
has no tractable unit cell (quasicrystals or amorphous met- 
als); they also facilitate dynamic or thermodynamic simula- 
tions (phonon spectra and phase transformations in complex 
alloys). 

A common modern approach to modelling atomic interac- 
tions classically is the "embedded-atom method" (EAM) [2] 
as well as "modified EAM" (3|], in which the full Hamilto- 
nian contains the usual pair term Vij (R), but also an implicitly 
many-atom term U(p), where p is a sum of contributions from 
nearby atoms. Accurate EAM potentials are straightforward 
to extract for monatomic systems but demand patience and 
skill to obtain even for binary systems 01; obviously, di- 
mensionality of parameter space becomes critical for multi- 
component systems. 

In this paper we report on an alternative approach [7] fitting 
only pair interactions but incorporating Friedel oscillations, 
optimized (typically) for a particular composition range. The 
oscillating analytical form is natural only for simple metals 
(Al, Mg,...) yet (see below) it works surprisingly well even 
when angular or many-body interactions are important. We 
will first describe the form of the potential and our methods for 
fitting it, then demonstrate its capabilities through case stud- 
ies in three alloy systems: Fe-B, Al-Mg-Zn, and Al-Cu-Fe; 
finally, we will summarize other systems where this method 
has been applied and discuss its limitations. 

Method: potential, database, and fitting — Our "empirical 
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All six parameters, including fc*, are taken as independent 
in the fit for each pair of elements. Eq. (Q]) was inspired by 
effective potentials (e.g. @] and Jgtl) used in previous work 
on structurally complex metals, e.g. quasicrystals ||7, 111. 
In such systems, energy differences between competing struc- 
tures are often controlled by second- and third-neighbor wells 
due to Friedel oscillations, which are a consequence (math- 
ematically) of Fourier transforming the Fermi surface, or 
(physically) are equivalent to the Hume-Rothery stabilization 
by enhancing the strength of structure factors that hybridize 
states across the Fermi surface. Ill 211 In the framework of 
Ref.[l2| (Sec. 6.6), the short-range repulsion is captured by the 
first term of (Q]); the medium-range potential (first-neighbor 
well) as well as the long-range oscillatory tail are captured by 
the second term, their relative weights being adjusted by the 
rji parameters. Note that empirically fitted potentials account 
for some of of the many-body contributions within their pair 
terms, which works better practically than truncating a sys- 
tematic expansion (e.g. GPT 11310 after the pair terms. 

The parameters in (Q]) are fitted to an ab-initio dataset (we 
always used the VASP code (lfo) combining both relaxed 
T = structures and molecular dynamics (MD) simula- 
tions at high T (usually the same structures, below the melt- 
ing point.) Structures are selected for the database within, or 
bracketing, the composition range of interest; when possible, 
a mix of simple and complex structures is used. A key crite- 
rion in choosing structures is to ensure an adequate number 
of contacts of each kind (in particular, nearest-neighbor con- 
tacts between the least abundant species). Also, all structures 
in the database should have similar atom densities fl411 . In 
particular, our high-T MD samples were constrained to have 
the same density as at T = 0, rather than the physical zero- 
pressure values. In MD simulations of the simpler structures, 
a supercell is always used with dimensions comparable to the 
potential cutoff radius, which (in this paper) is always 12 A. 

We define each structure's energy as a difference relative to 
a coexisting mixture (with the same total composition) of ref- 
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system 


B-Fe 


Al-Mg-Zn 


Al-Cu-Fe 


ref. 


Fe, BFe 2 


Al, Mg, Zn 


Al 7 Cu 2 Fe 








Al 2 Cu, Al 2 Fe 


data- 


BFe 2 [Al 2 Cu] 


MgZn 2 [Laves] 


Al 6 Fe 


base 


BFe 3 [CFe 3 ] 


T(Al 96 Mg 64 ) 


Al 5 Fe 2 [Al 23 CuFe 4 ] 




BFe 3 [Ni 3 P] 


[T(AlMgZn)] 


Al 2 Cu 




B 6 Fe 23 [C 6 Cr 23 ] 


T(Zn 9 6Mg 6 4) 


Al 7 Cu 2 Fe 




B 3 Fe 4 [B 3 Ni 4 ] 


[T(AlMgZn)] 


Al 23 CuFe 4 




BiFei [BCr] 


Ali 2 Mgi 7 


Al 72 Cu 4 Fe 24 [Al 3 Fe] 




a-Fe8oB 2 o 


AlMg4Znn 


"1/1-expe" 






/3'-Al 3 Mg 2 


"1/1-6D" 


Tm d 


1500 


450 (Al-Mg); 


500 ("1/1-6D"); 






1000 (Mg-Zn) 


1200 ("1/1-6D"); 








1000 (Al 7 Cu 2 Fe) 



TABLE I: Compounds used in each database used to fit potentials. 
Energies are expressed as difference from the tie-line or tie-plane 
defined by the simple reference structures listed; these structures are 
also included in the database. In brackets are the conventional names 
for crystal structure type. The temperature (in K) used for molecular- 
dynamics samples is Tmd (where compositions are given, those are 
the only T > inputs) 



erence phases, chosen to bracket all database compositions. 
Every structure is used for both forces (from MD at high T) 
and energy differences lfl5ll (high-T MD, as well as relaxed 
at T = 0). For the high-T portion, we took one snapshot 
of each structure at the end of a short ab-initio MD run. Typi- 
cally ~ 10 3 force components entered the fit, along with ~ 50 
energy differences (more for Al-Mg-Zn, fewer for Fe-B), and 
the forces are ~ 2-4 eV/A while the energy differences are 
~ 0.2-0.4 eV/atom; the fit residuals are ~ 5% and ~ 1% 
respectively. 

Our least-squares fit minimizes (by the Levenberg- 
Marquardt algorithm) x 2 = AEf/a E 2 + \ AF j? /^f 2 , 
where {A_E, } and {AF 3 } are the energy and force residuals; 
we found a weighting ratio ueI^f ~ 10~ 3 A was optimal 
so that neither energies nor forces dominate the fit. There is 
some risk of converging to a false minimum (or not converg- 
ing at all, from an unreasonable initial guess). Thus, it is im- 
portant to repeat the fit from several starting guesses. For this 
we used, e.g., potentials first fitted to pure elements or binary 
systems, and also used a library of parameter sets previously 
fitted for some different alloy system. The fitted parameters 
in (Q~|i for each of our examples are gathered in the table (TTJ 
similar potentials were plotted in Refs. [17] (for Al-Mg) and 
ifTsh (for Sc-Zn). 

Example 1: B-Fe. — Amorphous B-Fe is the simplest 
member of a family of technologically important metallic 
glasses. Our B-Fe database (Table included several crystals 
plus an "amorphous" sample with 100 atoms in an approxi- 
mately cubic box. The ab-initio data were calculated with a 
spin polarization so as to include magnetic contributions to 
the forces/energies. Il9ll For comparison, we also fitted the 
database to Lennard- Jones (LJ) potentials which lack oscilla- 
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FIG. 1: Scatter plot of pair potential result (vertical axis) versus ab- 
initio data (horizontal axis) for energy differences AEi (top panels) 
and forces Fj (bottom panels) The potentials are: (a,d) Lennard- 
Jones fit to B-Fe, (b,e) EOPP fit [Eq. {j}] to B-Fe, and (c,f) EOPP fit 
to Al-Mg-Zn. 



tions (starting this fit from the Kob-Anderson potentials [20]). 

We employed a modified tempering scheme (2lj] with 
replica exchange. That means two simulations were carried 
out for 100-atom samples of FesoB2o at T=1500K, one us- 
ing the ab-initio energy as its Hamiltonian and the other using 
EOPP. At intervals a Monte Carlo (MC) attempt is made to 
swap the two samples. These swaps accelerate the configura- 
tion sampling. The acceptance probability (we had ~ 70% for 
Fe-B) is a stringent test of how well the pair-potential Hamil- 
tonian mimics the ab-initio one. 

The superiority of the EOPP form over LJ is evident from 
the excellent diagnostic of scatter plots in Fig. Q] As an- 
other diagnostic, Fig. [2] shows the radial distribution functions 
found in an MD run at T =1500K, using the ab-initio ener- 
gies as well as both kinds of pair potential; the EOPP form 
faithfully reproduced all features of the ab-initio data. 

Example 2: Al-Mg-Zn — In this example, EOP poten- 
tials resolve site occupancies and the phase diagram of 
low-temperature modifications of the well-known cubic 



"Bergman" compound T(AlMgZn) 11221, 12311 (X-rays can 
hardly distinguish Al from Mg ions, as they differ by only 
one electron.) Our database (TablelJ) spanned the entire range 
of ternary compositions. (Al-Mg-Zn was previously modeled 
using pseudo-binary pair potentials, lHUl . - ) 

The EOP potentials were then used in a lattice-gas type MC 



annealing (using sites from diffraction II22L 12311 1 for each of 
91 compositions, with 52-64 Mg atoms and 24-72 Al atoms 
per cell (thus spanning the full Al-Zn range, and ±5% in Mg 
concentration). For 21 selected candidate structures, we sub- 
sequently did a full ab-initio calculation. Four of these struc- 
tures (Table Ullb were stable at T = 0. The first of these 
is the standard structure M22I1 (but with empty-centered icosa- 
hedral clusters); the second one differs by converting the Mg 
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FIG. 2: a,b,c: B-Fe pair potentials (c,d,e): Corresponding radial dis- 
tribution functions at 1500K, using ab-initio energies (dotted), em- 
pirical oscillating pair-potentials (solid), or fitted Lennard- Jones po- 
tentials (dashed). 



pair 


Cx 


m 


c 2 


7)2 k* if, 


B-B 


11.69 


5.184 


8.338 


5.149 3.950 2.617 


B-Fe 


2372.51 


17.232 


-6.043 


4.441 3.994 3.637 


Fe-Fe 


1.443 x10 s 


18.556 


4.939 


3.934 4.074 4.651 


Al-Al 


2669 


10.37 


-4.93 


4.41 3.88 -4.38 


Al-Mg 


4296 


10.61 


-39.4 


6.21 3.33 -4.02 


Mg-Mg 


14.59 


4.79 


183.3 


6.49 2.19 -4.53 


Al-Zn 


949 


9.36 


-2.52 


3.36 3.33 -2.16 


Mg-Zn 


282 


8.26 


-31.7 


5.58 2.82 -2.26 


Zn-Zn 


706 


9.33 


-6.01 


4.33 3.69 -3.48 


Al-Al 


4.431x10 s 


15.990 


-0.577 3.410 4.680 5.372 


Al-Fe 


1.842x10 s 


18.713 


7.884 


3.609 3.079 1.730 


Al-Cu 


3096.1 


12.079 


-1.457 3.261 3.687 2.978 


Fe-Fe 


125.26 


8.081 


1.115 


2.504 4.477 1.708 


Fe-Cu 


23.13 


5.856 


0.323 


1.520 3.086 1.947 


Cu-Cu 


11304 


12.631 


-2.124 3.301 2.830 0.159 



TABLE II: Parameters for B-Fe, Al-Mg-Zn and Al-Cu-Fe poten- 
tials [Eq. Q], with r in units of A and V(r) in eV. 



atoms with coordination 14 (the smallest for Mg) to Al atoms. 
We applied the same method to the recently solved /?'(AlMg) 
structure 11711 . a low-T variant of Samson's very large-cell 
/3(ATMg) iH; at T = /3'(AlMg) is essentially a line com- 
pound. 

Example 3: Al-Cu-Fe — The structures of the best thermo- 
dynamically stable icosahedral quasicrystals, e.g. i-Al-Mn- 
Pd and i-Al-Cu-Fe, are poorly known despite excellent long- 
range order and 20 years of study. This extends to a-AlCuFe, 
the so-called 1/1 cubic approximant to the quasicrystal (it is 
related to i-AlCuFe in the same way that T(AlMgZn) [see 
above] is related to i-AlMgZn.) 

Diffraction-based modeling, by itself, rarely resolves all the 
important structural details; this is apparent even in the a- 



composition 


Al(2) Mg(3) Zn(l) 


sp. gr. 


A_BvASP AiJpair 


Al 2 4Mg64Zn 7 2 


1 








Im3 


-122.1 


-117.3 


Al3eMg52Zn72 


1 


1 





Im3 


-114.1 


-109.0 


Al 4 8Mg 6 4Zn 4 8 


2/3 





1/3 


Immm 


-92.8 


-92.2 


Al 6 6Mg5gZn 3 6 


1/4 


1/2 





R3 


-77.9 


-76.1 



TABLE III: The four stable low-temperature modifications of the 
Bergman nhase TYAlMpZnl: the occunation is piven (as a fraction^ 
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FIG. 3: Partial phonon densities of states for Al, Cu and Fe:"calc" 
calculated Ck:alc") using EOP potentials, experimental ("exp") data 
after Ref. J3lh . 



AlCuFe crystal, for which the "solved" structure l25ll con- 
tains mixed-occupancy Al-Cu and Al-Fe sites; that uncer- 
tainty, is a substantial obstacle to realistic modelling of qua- 
sicrystal properties: the ab-initio energy of such a model is, 
we expect, unstable by ~400 meV/atom with respect to com- 
peting phases lf26ll ; the idealized and improved model 1 27] is, 
we found, unstable by 5 1 me V/atom. The root of the difficulty 
is that strong local correlations have not been accounted for, 
due to the averaging inherent in using Bragg peak intensities. 

We set up a robust database with 15312 force and 120 
energy datapoints (see Table H]). The fit converged to 0.2 
eV/A r.m.s. for forces and 8 me V/atom for energies. To 
refine the atomic structure of a(AlCuFe), we took can- 
didate sites from two diffraction-based models, either from 
a(AlCuFeSi) Q, or 1/1-AlCuRu and performed a 

lattice-gas Monte Carlo annealing using EOP potentials. The 
best configurations were then annealed by MD (still using 
EOPP) at moderately high temperatures (T =1000K) and fi- 
nally quenched. When their energies were recalculated ab- 
initio, the best example was still unstable, but only by the rel- 
atively small energy 35 me V/atom. The lowest energy struc- 
ture obtained as described above was then used for phonon 
DOS calculation by diagonalizing the dynamical matrix. The 
agreement between calculated and experimental DOS (Fig. 0) 
is encouraging. 

EOP potentials were also applied to predict the atomic 
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structure of the quasicrystal i-Al67.5Cu25Fei2.5, by simi- 
lar lattice-gas approach using (a periodic approximant of) a 
"quasilattice" of candidate sites. fl29ll The resulting structures 
were well-described by the formalism of a cut through a six- 
dimensional hyper-structure, but the hyper-atom occupancies 
differ in detail from the well-known Katz-Gratias model ll30ll . 
The ab-initio energy after relaxation was unstable (with re- 
spect to the competing phases [26]) by 45 meV/atom. The 
computed vibrational DOS agrees with neutron results 113 ill - 

Conclusion. — We have shown that empirical oscillating 
potentials with the simple analytic form (HJ mimic the atomic 
interactions of many metallic systems with sufficient accuracy 
to stand in for ab-initio energies when those would be compu- 
tationally prohibitive, e.g. to resolve mixed or partially occu- 
pied sites in diffraction refinements of complex structures, to 
measure thermodynamic expectations, or to compute phonon 
spectra. 

The easily obtained EOP potentials may serve as a starting 
point for the (tedious) construction of (more elaborate) EAM 
versions of the potentials. Furthermore, EOP potentials can 
accelerate ab-initio calculations, to identify uncompetitive al- 
ternatives to be ignored (as in our Al-Zn-Mg example) or to 
carry out all but the final iterations of a relaxation. EOP 
potentials appear to be fitted more robustly than EAM poten- 
tials for 4-component systems; we fitted B-C-Fe-Mo which 
enables for the first time MD simulations of technologically 
relevant metallic glasses (to probe the glass transition, two- 
level systems, etc.) . It should be noted, however, that any pair 
potential will fail in structures where the electron density has 
large variations in space (e.g. vacancies, edge dislocations, or 
surfaces). 

Rather than emphasize the easy cases (alloys rich in simple 
metals, e.g. AlMgZn, Mg-X, but also Zn-Y), we presented 
two difficult examples (B-Fe and Al-Cu-Fe), which are bor- 
derline for pair potentials, due to covalent bonding of nearest 
neighbor transition metals. In the many cases where bond di- 
rectionality is even more important, EOP potentials cannot be 
used. 

Our method had mixed success for pure elements: Al-Al 
potentials show an excellent fit to the ab-initio data, with Mg- 
Mg even better. But with pure Zn, pair potentials never re- 
produced the ab-initio phonon vibrational DOS; and the EOPP 
approach fails for pure Ga. 

EOP potentials show promise for many other alloy systems. 
They allowed extraction of long-wavelength phonons from 
molecular dynamics simulations IB3I1 in liquid Bi-Li system. 
Outstandingly accurate EOP were fitted for the system Mg- 
T (T=Pd or Ag) in the Mg-rich limit; for T=Ag, the Ag-Ag 
potential exhibits strong oscillations to large r that could only 
be fitted after the cutoff was extended to to 14A [32]. These 
systems include various complex phases based on icosahedral 
"Mackay" clusters (Mg42Xi2). EOP potentials have predicted 
- correctly - which Mg sites get substituted by Pd as the 
composition is varied 13211 in the phases 7(Mgo.sPdo.2) and 
5(Mgo.sPdo.2) [both "1/1 approximants", similar in structure 
to a-AlMnSi]. 
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